unit FloatMatrix;
{ some functions to do matrix and vector math , uses rm.na by default
  with NaNs as na
  Matrices are indexed with [col,row]
  }

interface

type
  TFloatVector = array of double;
  TFloatMatrix = array of array of double;

  TIntVector = array of integer;
  TIntMatrix = array of array of integer;

  TLogVector = array of boolean;
  TLogMatrix = array of array of boolean;

  TStringVector = array of string;
  TStringMatrix = array of array of string;

  function Log2IntVector(vector:TLogVector):TIntVector;

  function VectorCopy(V:TLogVector):TLogVector;overload;
  function MatrixCopy(M:TLogMatrix):TLogMatrix;overload;
  function GetRow(M:TLogMatrix;arow:integer):TLogVector;overload;
  function GetCol(M:TLogMatrix;acol:integer):TLogVector;overload;
  function SelectElements(V:TLogVector;vector:TIntVector):TLogVector;overload;
  function SelectElements(V:TLogVector;vector:TLogVector):TLogVector;overload;
  function SelectRows(M:TLogMatrix;vector:TIntVector):TLogMatrix;overload;
  function SelectRows(M:TLogMatrix;vector:TLogVector):TLogMatrix;overload;
  function SelectCols(M:TLogMatrix;vector:TIntVector):TLogMatrix;overload;
  function SelectCols(M:TLogMatrix;vector:TLogVector):TLogMatrix;overload;

  function VectorCopy(V:TIntVector):TIntVector;overload;
  function MatrixCopy(M:TIntMatrix):TIntMatrix;overload;
  function GetRow(M:TIntMatrix;arow:integer):TIntVector;overload;
  function GetCol(M:TIntMatrix;acol:integer):TIntVector;overload;
  function SelectElements(V:TIntVector;vector:TIntVector):TIntVector;overload;
  function SelectElements(V:TIntVector;vector:TLogVector):TIntVector;overload;
  function SelectRows(M:TIntMatrix;vector:TIntVector):TIntMatrix;overload;
  function SelectRows(M:TIntMatrix;vector:TLogVector):TIntMatrix;overload;
  function SelectCols(M:TIntMatrix;vector:TIntVector):TIntMatrix;overload;
  function SelectCols(M:TIntMatrix;vector:TLogVector):TIntMatrix;overload;

  function rm_na(V:TFloatVector):TFloatVector;overload;
  function VectorCopy(V:TFloatVector):TFloatVector;overload;
  function MatrixCopy(M:TFloatMatrix):TFloatMatrix;overload;
  function GetRow(M:TFloatMatrix;arow:integer):TFloatVector;overload;
  function GetCol(M:TFloatMatrix;acol:integer):TFloatVector;overload;
  function SelectElements(V:TFloatVector;vector:TIntVector):TFloatVector;overload;
  function SelectElements(V:TFloatVector;vector:TLogVector):TFloatVector;overload;
  function SelectRows(M:TFloatMatrix;vector:TIntVector):TFloatMatrix;overload;
  function SelectRows(M:TFloatMatrix;vector:TLogVector):TFloatMatrix;overload;
  function SelectCols(M:TFloatMatrix;vector:TIntVector):TFloatMatrix;overload;
  function SelectCols(M:TFloatMatrix;vector:TLogVector):TFloatMatrix;overload;

  function VectorCopy(V:TStringVector):TStringVector;overload;
  function MatrixCopy(M:TStringMatrix):TStringMatrix;overload;
  function GetRow(M:TStringMatrix;arow:integer):TStringVector;overload;
  function GetCol(M:TStringMatrix;acol:integer):TStringVector;overload;
  function SelectElements(V:TStringVector;vector:TIntVector):TStringVector;overload;
  function SelectElements(V:TStringVector;vector:TLogVector):TStringVector;overload;
  function SelectRows(M:TStringMatrix;vector:TIntVector):TStringMatrix;overload;
  function SelectRows(M:TStringMatrix;vector:TLogVector):TStringMatrix;overload;
  function SelectCols(M:TStringMatrix;vector:TIntVector):TStringMatrix;overload;
  function SelectCols(M:TStringMatrix;vector:TLogVector):TStringMatrix;overload;

  function Mean(V:TFloatVector):double;overload;
  function Mean(M:TFloatMatrix):double;overload;
  function Mean(M:TFloatMatrix;margin:integer):TFloatVector;overload;

  function Median(V:TFloatVector):double;overload;
  function Median(M:TFloatMatrix):double;overload;
  function Median(M:TFloatMatrix;margin:integer):TFloatVector;overload;

  function Sort(v:TFloatVector):TFloatVector;
  function SortIndex(v:TFloatVector):TIntVector;

  function Seq(StartValue,EndValue:integer):TIntVector;

implementation
uses IEEE754,Math;

function Log2IntVector(vector:TLogVector):TIntVector;
var
  i:integer;
begin
  Setlength(result,0);
  for i:=0 to High(vector) do
    if vector[i] then
    begin
      SetLength(result,Length(result)+1);
      result[High(result)]:=i;
    end;
end;

function Matrix2Vector(M:TFloatmatrix):TFloatVector;
var
  i,j : integer;
begin
  SetLength(result,length(M)*length(M[0]));
  for i:=0 to High(M) do
    for j:=0 to High(M[0]) do
      result[i * length(M) + j]:=M[i,j];
end;

function rm_na(V:TFloatVector):TFloatVector;
var
  i,valid:integer;
begin
  valid:=Low(V);
  for i:=Low(V) to High(V) do
    if Not IsNan(V[i]) then
    begin
      V[valid]:=V[i];
      inc(valid);
    end;
  SetLength(v,valid);
  result:=v;
end;

{ Log Subsetting and copy}
function VectorCopy(V:TLogVector):TLogVector;overload;
var
  elt:integer;
begin
  Setlength(result,Length(V));
  for elt:=Low(V) to High(V) do
    result[elt]:=V[elt];
end;

function MatrixCopy(M:TLogMatrix):TLogMatrix;
var
  acol,arow:integer;
begin
  Setlength(result,Length(M),Length(M[0]));
  for acol:=Low(M) to High(M) do
    for arow:=Low(M[0]) to High(M[0]) do
      result[acol,arow]:=M[acol,arow];
end;

function GetRow(M:TLogMatrix;arow:integer):TLogVector;
var
  acol:integer;
begin
  SetLength(result,Length(M));
  for acol:=Low(M) to High(M) do
    result[acol]:=M[acol,arow];
end;

function GetCol(M:TLogMatrix;acol:integer):TLogVector;
var
  arow:integer;
begin
  SetLength(result,Length(M[acol]));
  for arow:=Low(M[acol]) to High(M[acol]) do
    result[arow]:=M[acol,arow];
end;

function SelectElements(V:TLogVector;vector:TIntVector):TLogVector;overload
var
  elt:integer;
begin
  Setlength(result,Length(vector));
  for elt:=0 to High(Vector) do
    result[elt]:=V[Vector[elt]];
end;

function SelectElements(V:TLogVector;vector:TLogVector):TLogVector;overload
begin
  if Length(vector)=Length(V) then
    result:=SelectElements(V,Log2Intvector(Vector))
  else
    raise EInvalidArgument.Create('Length of logical vector does not match number of elements in vector');
end;

function SelectRows(M:TLogMatrix;vector:TIntVector):TLogMatrix;overload;
var
  acol,arow:integer;
begin
  Setlength(result,Length(M),Length(vector));
  for acol:=0 to High(M) do
    for arow:=0 to High(Vector) do
      result[acol,arow]:=M[acol,Vector[arow]];
end;

function SelectRows(M:TLogMatrix;vector:TLogVector):TLogMatrix;overload;
begin
  if Length(vector)=Length(M[0]) then
    result:=SelectRows(M,Log2Intvector(Vector))
  else
    raise EInvalidArgument.Create('Length of logical vector does not match number of rows in matrix');
end;

function SelectCols(M:TLogMatrix;vector:TIntVector):TLogMatrix;overload;
var
  acol,arow:integer;
begin
  Setlength(result,Length(vector),Length(M[0]));
  for acol:=0 to High(Vector) do
    for arow:=0 to High(M[0]) do
      result[acol,arow]:=M[Vector[acol],arow];
end;

function SelectCols(M:TLogMatrix;vector:TLogVector):TLogMatrix;overload;
begin
  if Length(vector)=Length(M) then
    result:=SelectCols(M,Log2Intvector(Vector))
  else
    raise EInvalidArgument.Create('Length of logical vector does not match number of columns in matrix');
end;

{ Int Subsetting and copy}
function VectorCopy(V:TIntVector):TIntVector;overload;
var
  elt:integer;
begin
  Setlength(result,Length(V));
  for elt:=Low(V) to High(V) do
    result[elt]:=V[elt];
end;

function MatrixCopy(M:TIntMatrix):TIntMatrix;
var
  acol,arow:integer;
begin
  Setlength(result,Length(M),Length(M[0]));
  for acol:=Low(M) to High(M) do
    for arow:=Low(M[0]) to High(M[0]) do
      result[acol,arow]:=M[acol,arow];
end;

function GetRow(M:TIntMatrix;arow:integer):TIntVector;
var
  acol:integer;
begin
  SetLength(result,Length(M));
  for acol:=Low(M) to High(M) do
    result[acol]:=M[acol,arow];
end;

function GetCol(M:TIntMatrix;acol:integer):TIntVector;
var
  arow:integer;
begin
  SetLength(result,Length(M[acol]));
  for arow:=Low(M[acol]) to High(M[acol]) do
    result[arow]:=M[acol,arow];
end;

function SelectElements(V:TIntVector;vector:TIntVector):TIntVector;overload
var
  elt:integer;
begin
  Setlength(result,Length(vector));
  for elt:=0 to High(Vector) do
    result[elt]:=V[Vector[elt]];
end;

function SelectElements(V:TIntVector;vector:TLogVector):TIntVector;overload
begin
  if Length(vector)=Length(V) then
    result:=SelectElements(V,Log2Intvector(Vector))
  else
    raise EInvalidArgument.Create('Length of logical vector does not match number of elements in vector');
end;

function SelectRows(M:TIntMatrix;vector:TIntVector):TIntMatrix;overload;
var
  acol,arow:integer;
begin
  Setlength(result,Length(M),Length(vector));
  for acol:=0 to High(M) do
    for arow:=0 to High(Vector) do
      result[acol,arow]:=M[acol,Vector[arow]];
end;

function SelectRows(M:TIntMatrix;vector:TLogVector):TIntMatrix;overload;
begin
  if Length(vector)=Length(M[0]) then
    result:=SelectRows(M,Log2Intvector(Vector))
  else
    raise EInvalidArgument.Create('Length of logical vector does not match number of rows in matrix');
end;

function SelectCols(M:TIntMatrix;vector:TIntVector):TIntMatrix;overload;
var
  acol,arow:integer;
begin
  Setlength(result,Length(vector),Length(M[0]));
  for acol:=0 to High(Vector) do
    for arow:=0 to High(M[0]) do
      result[acol,arow]:=M[Vector[acol],arow];
end;

function SelectCols(M:TIntMatrix;vector:TLogVector):TIntMatrix;overload;
begin
  if Length(vector)=Length(M) then
    result:=SelectCols(M,Log2Intvector(Vector))
  else
    raise EInvalidArgument.Create('Length of logical vector does not match number of columns in matrix');
end;

{ Float Subsetting and copy}
function VectorCopy(V:TFloatVector):TFloatVector;overload;
var
  elt:integer;
begin
  Setlength(result,Length(V));
  for elt:=Low(V) to High(V) do
    result[elt]:=V[elt];
end;

function MatrixCopy(M:TFloatMatrix):TFloatMatrix;
var
  acol,arow:integer;
begin
  Setlength(result,Length(M),Length(M[0]));
  for acol:=Low(M) to High(M) do
    for arow:=Low(M[0]) to High(M[0]) do
      result[acol,arow]:=M[acol,arow];
end;

function GetRow(M:TFloatMatrix;arow:integer):TFloatVector;
var
  acol:integer;
begin
  SetLength(result,Length(M));
  for acol:=Low(M) to High(M) do
    result[acol]:=M[acol,arow];
end;

function GetCol(M:TFloatMatrix;acol:integer):TFloatVector;
var
  arow:integer;
begin
  SetLength(result,Length(M[acol]));
  for arow:=Low(M[acol]) to High(M[acol]) do
    result[arow]:=M[acol,arow];
end;

function SelectElements(V:TFloatVector;vector:TIntVector):TFloatVector;overload
var
  elt:integer;
begin
  Setlength(result,Length(vector));
  for elt:=0 to High(Vector) do
    result[elt]:=V[Vector[elt]];
end;

function SelectElements(V:TFloatVector;vector:TLogVector):TFloatVector;overload
begin
  if Length(vector)=Length(V) then
    result:=SelectElements(V,Log2Intvector(Vector))
  else
    raise EInvalidArgument.Create('Length of logical vector does not match number of elements in vector');
end;

function SelectRows(M:TFloatMatrix;vector:TIntVector):TFloatMatrix;overload;
var
  acol,arow:integer;
begin
  Setlength(result,Length(M),Length(vector));
  for acol:=0 to High(M) do
    for arow:=0 to High(Vector) do
      result[acol,arow]:=M[acol,Vector[arow]];
end;

function SelectRows(M:TFloatMatrix;vector:TLogVector):TFloatMatrix;overload;
begin
  if Length(vector)=Length(M[0]) then
    result:=SelectRows(M,Log2Intvector(Vector))
  else
    raise EInvalidArgument.Create('Length of logical vector does not match number of rows in matrix');
end;

function SelectCols(M:TFloatMatrix;vector:TIntVector):TFloatMatrix;overload;
var
  acol,arow:integer;
begin
  Setlength(result,Length(vector),Length(M[0]));
  for acol:=0 to High(Vector) do
    for arow:=0 to High(M[0]) do
      result[acol,arow]:=M[Vector[acol],arow];
end;

function SelectCols(M:TFloatMatrix;vector:TLogVector):TFloatMatrix;overload;
begin
  if Length(vector)=Length(M) then
    result:=SelectCols(M,Log2Intvector(Vector))
  else
    raise EInvalidArgument.Create('Length of logical vector does not match number of columns in matrix');
end;

{ String Subsetting and copy}
function VectorCopy(V:TStringVector):TStringVector;overload;
var
  elt:integer;
begin
  Setlength(result,Length(V));
  for elt:=Low(V) to High(V) do
    result[elt]:=V[elt];
end;

function MatrixCopy(M:TStringMatrix):TStringMatrix;
var
  acol,arow:integer;
begin
  Setlength(result,Length(M),Length(M[0]));
  for acol:=Low(M) to High(M) do
    for arow:=Low(M[0]) to High(M[0]) do
      result[acol,arow]:=M[acol,arow];
end;

function GetRow(M:TStringMatrix;arow:integer):TStringVector;
var
  acol:integer;
begin
  SetLength(result,Length(M));
  for acol:=Low(M) to High(M) do
    result[acol]:=M[acol,arow];
end;

function GetCol(M:TStringMatrix;acol:integer):TStringVector;
var
  arow:integer;
begin
  SetLength(result,Length(M[acol]));
  for arow:=Low(M[acol]) to High(M[acol]) do
    result[arow]:=M[acol,arow];
end;

function SelectElements(V:TStringVector;vector:TIntVector):TStringVector;overload
var
  elt:integer;
begin
  Setlength(result,Length(vector));
  for elt:=0 to High(Vector) do
    result[elt]:=V[Vector[elt]];
end;

function SelectElements(V:TStringVector;vector:TLogVector):TStringVector;overload
begin
  if Length(vector)=Length(V) then
    result:=SelectElements(V,Log2Intvector(Vector))
  else
    raise EInvalidArgument.Create('Length of logical vector does not match number of elements in vector');
end;

function SelectRows(M:TStringMatrix;vector:TIntVector):TStringMatrix;overload;
var
  acol,arow:integer;
begin
  Setlength(result,Length(M),Length(vector));
  for acol:=0 to High(M) do
    for arow:=0 to High(Vector) do
      result[acol,arow]:=M[acol,Vector[arow]];
end;

function SelectRows(M:TStringMatrix;vector:TLogVector):TStringMatrix;overload;
begin
  if Length(vector)=Length(M[0]) then
    result:=SelectRows(M,Log2Intvector(Vector))
  else
    raise EInvalidArgument.Create('Length of logical vector does not match number of rows in matrix');
end;

function SelectCols(M:TStringMatrix;vector:TIntVector):TStringMatrix;overload;
var
  acol,arow:integer;
begin
  Setlength(result,Length(vector),Length(M[0]));
  for acol:=0 to High(Vector) do
    for arow:=0 to High(M[0]) do
      result[acol,arow]:=M[Vector[acol],arow];
end;

function SelectCols(M:TStringMatrix;vector:TLogVector):TStringMatrix;overload;
begin
  if Length(vector)=Length(M) then
    result:=SelectCols(M,Log2Intvector(Vector))
  else
    raise EInvalidArgument.Create('Length of logical vector does not match number of columns in matrix');
end;

function Mean(V:TFloatVector):double;overload;
var
  i:integer;
begin
  V:=rm_na(V);
  if Length(V)>0 then
  begin
    for i:=1 to High(V) do
      V[0]:=V[0]+V[i];
    result:=V[0] / Length(V);
  end
  else
    result:=NaN;
end;

function Mean(M:TFloatMatrix):double;overload;
begin
  Result:=Mean(Matrix2Vector(M));
end;

function Mean(M:TFloatMatrix;margin:integer):TFloatVector;overload;
var
  i:integer;
begin
  case margin of
  1:begin
    Setlength(result,Length(M));
    for i:=Low(M) to High(M) do
      result[i]:=Mean(GetCol(M,i));
  end;
  2:begin
    Setlength(result,Length(M[0]));
    for i:=Low(M[0]) to High(M[0]) do
      result[i]:=Mean(GetRow(M,i));
  end;
  end;
end;


function Median(V:TFloatVector):double;overload;
var
  n:integer;
begin
  V:=rm_na(V);
  V:=Sort(V);
  n:=length(V);
  if n>0 then
  begin
    if Odd(n) then
      result:=V[n div 2]
    else
      result:=(V[(n div 2) -1] + V[n div 2]) / 2;
  end
  else
    result:=NaN;
end;

function Median(M:TFloatMatrix):double;overload;
begin
  result:=Median(Matrix2Vector(M));
end;

function Median(M:TFloatMatrix;margin:integer):TFloatVector;overload;
var
  i:integer;
begin
  case margin of
  1:begin
    Setlength(result,Length(M));
    for i:=Low(M) to High(M) do
      result[i]:=Median(GetCol(M,i));
  end;
  2:begin
    Setlength(result,Length(M[0]));
    for i:=Low(M[0]) to High(M[0]) do
      result[i]:=Median(GetRow(M,i));
  end;
  end;
end;


function Sort(v:TFloatVector):TFloatVector;
var
  i,j,n : integer;
  tmp:double;

begin
  n:=length(v);
  for i:=0 to n-2 do
    for j:=i+1 to n-1 do
      if Not isNan(v[j]) and (IsNan(v[i]) or (v[i]>v[j])) then
      begin
        tmp:=v[j];
        v[j]:=v[i];
        v[i]:=tmp;
      end;
  result:=v;
end;

function SortIndex(v:TFloatVector):TIntVector;
var
  i,j,n : integer;
  tmpval:double;
  tmpidx:integer;
begin
  n:=length(v);
  Setlength(result,n);
  for i:=0 to n-1 do
    result[i]:=i;
  for i:=0 to n-2 do
    for j:=i+1 to n-1 do
      if Not isNan(v[j]) and (IsNan(v[i]) or (v[i]>v[j])) then
      begin
        tmpval:=v[j];
        v[j]:=v[i];
        v[i]:=tmpval;
        tmpidx:=result[j];
        result[j]:=result[i];
        result[i]:=tmpidx;
      end;
end;

function Seq(StartValue,EndValue:integer):TIntVector;
var
  by,i:integer;
begin
  SetLength(result,abs(EndValue-StartValue)+1);
  if StartValue>EndValue then
    by:=-1
  else
    by:=1;
  Endvalue:=Abs(EndValue-StartValue);
  for i:=0 to Endvalue do
  begin
    result[i]:=StartValue;
    Inc(StartValue,by);
  end;
end;

end.
